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ABSTRACT 

The basic effect of the earth's atmospheric refraction on telescope operation is the reduction 
of the true zenith angle to the apparent zenith angle, associated with prismatic aberrations 
due to the dispersion in air. If one attempts coherent superposition of star images in ground- 
based interferometry, one is in addition interested in the optical path length associated with 
the refracted rays. In a model of a flat earth, the optical path difference between these is not 
concerned as the translational symmetry of the setup means no net effect remains. 

Here, I evaluate these interferometric integrals in the more realistic arrangement of two tele- 
scopes located on the surface of a common earth sphere and point to a star through an atmosphere 
which also possesses spherical symmetry. Some focus is put on working out series expansions in 
terms of the small ratio of the baseline over the earth radius, which allows to bypass some nu- 
merics which otherwise is challenged by strong cancellation effects in building the optical path 
difference. 

Subject headings: Atmospheric Refraction; Optical Interferometry; Baseline; Spherical Geometry; Zenith 



1. Pointing and Apparent Star Altitudes 
1.1. Basics: Flat Earth Model 

The standard model of an interferometric setup 
and delay line correction for a star at the true 
zenith angle z is shown in Fig. 1: The optical path 
difference (OPD) D shows up once in the vacuum 
above telescope 1, and is added for telescope 2 
on the ground at some local index of refraction. 
The atmosphere is horizontally homogeneous and 
the earth flat; therefore no correction is needed 
for the ray path curvature induced by any verti- 
cal gradient of the index of refraction through the 
atmosphere, because these two paths match each 
other at each height above ground. 

The pointing difference R between the appar- 
ent and actual altitude of a star that is induced by 
the refraction of the earth atmosphere is in a sim- 
ple model of a flat earth (Green 1985; Filippcnko 
1982) 

R ss (no — 1) tan zq, (1) 



where no is the index of refraction on the ground, 
zq the apparent zenith distance on the ground, and 



R = z - z > 



(2) 



obtained in radian. 

This implies chromatic effects (Livengood et al. 
1999; Langhans et al. 2003; Colavita et al. 2004; 
Roe 2002) and dependencies on the atmospheric 
model, commonly summarized under the label 
"transversal atmospheric dispersion" in Astron- 
omy. There is a rainbow effect in Eq. (1): as no 
is a function of the wavelength, R becomes dis- 
persive, too: between wavelengths of 2 /im and 
2.4 /im, we get a difference of Ano « 1.05 • 10~ 7 
at 2600 m above sea level, which translates into a 
spectral smear of AR ss 22 mas • tanz . 

Eq. (1) suggests that the relative error in the 
star's altitude definition is close to the relative er- 
ror in the dielectric function and susceptibility at 
the telescope site. In addition, Snell's law of re- 
fraction states that product nsinz between the 
height-dependent index of refraction n and the 
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Fig. 1 . — The standard model of delay correction 
and recovery for a star at zenith angle z observed 
by two telescopes that look through horizontal lay- 
ers of the stratified atmosphere. D = frsinz = 
Ptanz. 

sine of the angle of refraction at that height re- 
mains constant along the path (Green 1985, (4.2)). 
Therefore the optical path delay is D = 6sinz = 
6nosinz ; it changes as the astronomical object 
changes position in z, or, supposed b is fixed, ac- 
cording to atmospheric parameters accessible at 
the ground level. The benefit of this analysis is 
that both, the pointing correction R in Eq. (1) and 
the OPD measured on the ground, are functions 
of the index of refraction at the telescope site, not 
functionals of the entire layered atmosphere. 

The theme of this paper are corrections to 
these statements considering a non-turbulent at- 
mosphere covering an earth surface of constant, 
but non-negligible curvature. The rest of Sec. 1 
shortly describes the standard theory of refraction 
and defines two different baseline lengths. Sec. 
2 concentrates on the integral formulation of the 
OPD calculation through the atmosphere: geome- 
tries with constrained azimuths suffice to intro- 
duce all relevant concepts; general star positions 
are then reduced to the constrained case. 

1.2. Spherical Earth: Geometry 

Things are more complicated if we start to look 
at the more realistic model of a spherical earth. A 
telescope distance of b — 100 m on an earth of ra- 



dius p = 6368 km leads to a pointing mismatch of 
purely geometric origin of about bj p w 3.2 arcsec 
(Fig. 2). The baseline b is the distance between 
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Fig. 2. — Two telescopes placed on the earth of 
radius p at a distance b looking at the same star 
experience local zenith angles that differ by Z rts 
b/p rad. 

the telescope locations on the earth, the length of 
the straight secant drawn in Figs. 2 and 3, 



b = /V 2 (! -cosZ) 
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The inversion of this series reads 

■bV _3_ 

s p) + 640 



b _1_ 

~p + 24 



(3) 
(4) 

(5) 



(6) 



The angle approximation Z w b/p is an estimate 
to b/p = 2sin(Z/2), a limit of a baseline so short 
that it does not matter whether it is measured 
along a straight line (as drawn in Fig. 2) or along 
the circular perimeter. The relative error in this 
approximation is « Z 2 /24, or w 10~ n for the ex- 
ample of b — 100 m. 

Pointing/guiding is a functionality of the indi- 
vidual telescopes: the existence of a nonzero point- 
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ing difference is absorbed in the telescope opera- 
tion, and any delay originating from there is to 
first order recovered in the tracking. We hereby 
explicitly discard any "separation" term of Gublcr 
and Tytler eminent from their consideration of 
single telescopes (Gubler & Tytler 1998), and ac- 
knowledge that both telescopes are "mispointing" 
at the same time. 

This cross-eyed geometry indirectly transforms 
into a contribution to the delay that may be 
understood — and calculated to lowest order — on 
the basis of what is said in Sec. 1.1, but has no par- 
allel with single telescopes as long as their diame- 
ters are much smaller than the earth radius: The 
example of 3.2 arcsec from above translates into 
an additional angle of refraction of AR w Azo(l + 
tan 2 z )(n - 1) > Az (n - 1) w Z(n - 1) w 650 
/xas at 2600 m above sea level, which is the first 
derivative of Eq. (1) w.r.t. z . Dropping the term 
tan 2 Zo here means this is a lower estimate in the 
limit of stars at the zenith. 

The effect on the delay could be understood 
in the "standard model of delay line correction," 
where the two rays of the star that will eventually 
hit the two telescopes "generate" a phase differ- 
ence in the vacuum (undone later on in the de- 
lay line tunnel) as they hit the top layer of the 
earth atmosphere with a path difference D (Fig. 
1). The curvature correction means that this top 
layer "bends back" a bit more for the telescope fur- 
ther away from the star, which slightly increases 
the angle of incidence on the atmosphere. 

Fig. 3 shows that the only obvious definition of 
the baseline length b is on earth. Three alterna- 
tives have been marked with question marks in the 
figure: 

1. Simple outward projection along some "com- 
mon" zenith direction of both telescopes is 
futile, because it would only scale b with 
some arbitrary distance from the earth cen- 
ter. 

2. Some sliding definition along the curved rays 
fails because there are no fix points or right 
angles for reference. 

3. The baseline cannot be rigorously defined far 
outside the atmosphere. Actually only the 
projection P of the two telescopes remains 
well defined as the rays can be considered 
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■ earth fixed 
coordinate system 

D 




Fig. 3. — Interferometry with two telescopes lo- 
cated on a non-planar surface needs to be aware 
of the vagueness of the definition of a baseline b. 

parallel, and so remains some direction angle 
a relative to some geostationary coordinate 
system. One could try to use the b? in Fig. 
3 as some line parallel to the earth-based b 
anchored at P, but needs to stay aware of 
the fact that this replaces the zenith angle 
z in Fig. (1) by the angle (5 = a — <fi, where 
a is the polar coordinate of the star and <f> 
the mean polar coordinate of the telescopes 
in the geocentric coordinate system. 

In a formal way, I define an effective baseline 
length b* above the atmosphere via 

b* 2 = P 2 + D 2 (7) 

which turns out to be longer than b because the 
earth's atmosphere represents a gradient lens (Hui 
& Seager 2002). This effective baseline is a func- 
tion of the baseline b and of the star zenith an- 
gle, and a functional of the atmospheric refraction 
n(r). 

1.3. Spherical Earth: Atmospheric Layers 

The delay line modification described in the 
previous paragraph is of purely geometrical nature 
and is closely related to the discussion of differ- 
ences between two light rays, but there are more 
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pointing "corrections" that already show up in the 
single ray case, as discussed next. 

Eq. (1) is based on the assumption that the gra- 
dient of the refractive index n along the ray path is 
parallel to a global, constant zenith vector. As the 
gradient and the zenith vector change in direction 
along the path through a spherically symmetric 
atmosphere, and as we assume that n becomes a 
function of the radial distance to the spherical sur- 
face of the earth, the equation becomes more ac- 
curately (Green 1985, (4. 19)) (Thomas & Joseph 
1995; Aucr & Standish 2000; Nener et al. 2003; 
Noerdlinger 1999; Tannous & Nigrin 2001) 



R = priQ sin zq 



i; 



dn 



n(r 2 n 2 — p 2 n\ sin 2 



(8) 

an integral over the refractive index, to start above 
the atmosphere (n = 1) and to end at the telescope 
position (n = no), and where r > p is the distance 
to the earth center. Its Taylor expansion is com- 
monly written as an expansion in powers of tan zo 
(Stone 1996), the observed quantity, 



R = 



r na dn 
pn tan z / — 



n*r 



+ 



+ 



pn tan 3 z f 

2 h 
3pno tan 5 zq 



no p 2 n 2 _ r 2 n 2 



dn 



i: 



(p 2 nl - r 2 n 2 ) 2 



-dn + .{9) 



The major new aspect is that the pointing direc- 
tion becomes a functional of the height spectrum 
of the index n, and the angle of arrival becomes 
site-dependent (Conan et al. 2000). 

Accurate modeling of n(r) is not within the 
scope of this treatise here and unambitiously dis- 
cussed in App. A. Tables and graphs to follow use 
an exponential depletion of the susceptibility 



X 



1 



1 



(10) 



to the vacuum of the universe with a scale height 
K, 

x(r)=Xoe- {r - p)/K . (11) 
The explicit parameters are 

Xo = 4 • 10~ 4 , p = 6380km, K = 10km, (12) 

unless otherwise noted, representing a prototyp- 
ical K-band value at 2600 m above sea level 



(Mathar 2004). Within this model, Fig. 4 shows 
the change in R introduced by switching from the 
flat earth model, 



R — > n sin z 



i: 



dn 



n(n 2 — n\ sin 2 zq) 1 ^ 2 ' 
to the spherical model of the atmosphere. 



(13) 
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Fig. 4. — Solid line: The relative error in the re- 
fraction angle R, in percent, made by switching 
from the spherical model, Eq. (8), to the limit of 
the flat earth, Eq. (13). Dashed line: The relative 
error by switching from the spherical model to the 
approximation (1). The absolute errors approach 
as z — > 0. 



2. Accumulated Optical Path Lengths for 
Spherical Geometry 

2.1. Planar, Overhead Geometry 

2.1.1. Single Star 

The optical path length along the curved ray 
trajectory is the line integral of the refractive in- 
dex over the geometric path, similar to Eq. (8) 



Jn=\ 



n ° dr 
n- 



cosi/> 



(14) 



where dr j cos ip is the length of the diagonal path 
element in Fig. 5. ip the local zenith angle of the 
star, and the r the distance from the earth cen- 
ter (Green 1985, Fig. 4.4). This quantity includes 
the geometric path length up to the star and is 
infinite in our applications. We assume a maxi- 
mum height H of the atmosphere, n\ r>p+H = 1, 
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Fig. 5. — The integrals (8) and (14) accumu- 
late the change in the geocentric zenith angle 
ip{r) and the optical path length n/ cosip, respec- 
tively, along the light path, applying Snell's law 
rn sin tp(r) = const at each differential layer. 



with the option to look at the limit H — > oo if the 
density has no clear ceiling, like in the models of 
Eq. (11). We define the impact parameters I\ and 
I 2 of the rays, < I2 < I\ < p, which were their 
smallest distances to the earth center if they would 
pass by along geometric straight lines without any 
diffraction — as used in atomic collision theory. 



P = h-h 



(15) 



is the projected baseline, measured above the at- 
mosphere. Sec. 2.1 considers the case in which the 
star, the two telescopes (the baseline), and the 
earth center are coplanar: Fig. 6. Comments on 
the general configuration of an unconstrained star 
azimuth follow in Sec. 2.2. 

The difference between the integrals Eq. (14) 



D = L x - L 2 



I 



r—p 



ndr 



P+H 



,dr 



cos ^( 2 ) 



(16) 



(17) 



r—p 



is the OPD for two telescopes at a common base- 
line b, with different apparent and different true 



Fig. 6. — Sketch of the impact parameters h, the 

(i) 

auxiliary geocentric zenith angles ip H ■> the true 
zenith angles Zi, the angles of refraction Ri, and 
the projected baseline P. The two light rays en- 
ter from the right, hit the upper atmosphere at a 
height H where indicated by the dotted quarter- 
circle, and eventually the telescopes Ti. 

zenith angles (individual pointing), but looking at 
the same star. The additional straight line seg- 
ment for the ray to telescope 1 through vacuum 
before it reaches the altitude H is 



D v = ^{p + HY-P 2 -^( fi + HY-ll 

= {p + H)\<xxi%l$ -cosVS 5 ]- (18) 

The Taylor series of this term up to third order in 
P is 



D v w Ptan^ ) + 



tanV^ ) _ p2+ _ 



(2) 



1 ta,n 6 ip y * 



2-^2 COS 2 2 J| CO g2 ^(2) 

(19) 

We separate this piece from the integral for the 
telescope 1, 

r p+H 



D = D V 



cosi/^ 1 ' cos?/^ 2 ) 



dr, (20) 



where the major difference w.r.t. Fig. 1 is that this 
integral does not vanish, because the atmosphere 
is now hit at two different angles ipjp 7^ ■ 
The integral would also not vanish, if the factor 
n would be dropped to deduce the geometric path 
difference of the curved beams. 

The Constance of the product rnsimp along 
each curved trajectory (Green 1985, (4.16)), 

rnsinV> (1) = (p + H) sin = I\ ; (21) 
rnsinV> (2) = (p + H) sin^ ) = h, (22) 
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is inserted into the previous equation, 

P+H 

D = D v + J rn 

r—p 



i 



yW 2 - {p + H) 2 sin 2 ^ 
dr{23) 



1 



r 2 n 2 -(p + H) 2 sin 2 ^ 



The term in square brackets allows another 
Taylor expansion 



1 



1 



y r 2 n 2 



:*n 2 - if 



p | I r 2 n 2 + 2I 2 ^ p2 



(r 2 n 2 — I 2 ) 3 / 2 2 (r 2 n 2 - I 2 ) 5 / 2 ' 



1 3r 2 n 2 
+ - 



211 



2(r 2 n 2 -I 2 )V- 



I 2 P 3 - 



(24) 



All correction terms of the spherical geometry in 
Eqs. (19)-(24) have a positive sign. The contri- 
bution of the term of 0(P 2 ) amounts to w 2 mm, 
and of the term of 0(P 3 ) to w 60 nm as the zenith 
angle approaches 60 deg at b = 100 m. A consis- 
tency check of Eq. (24) is that its contribution of 
the first, linear Taylor order to the integral in Eq. 
(23) becomes 



P+H 



hP 



I r—p 
n->l 



(r 2 n 2 - I 2 ) 3 1 2 



dr 



Pt&niP {2) +P 



h 



T 2 



(25) 



in the vacuum limit, such that the term Ptanipff 
cancels the first term in Eq. (19), and only the 
term Ptanz depending on the true topocentric 
zenith angle remains, equivalent to Fig. 1. 

The influence of the spherical geometry on 
the atmospheric path length difference is demon- 
strated in Fig. 7. 



2.1A 



Astrometry (two true zenith distances) 



Up to here, an accurate computation of the de- 
lay D by integration over the atmosphere layers 
is not competitive against measuring the equiva- 
lent value on the ground (Fig. 1), since the gas 
densities on the ground along the beam path — 
input to calculation of the refractive indexes — are 
accessible to sensors and much better known than 
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Fig. 7. — The difference between the full integral 
(23) of the delay D and the approximation D w 
bsinz, which one would derive from Fig. 1 using 
the mean true zenith angle z = (z\ + Z2)/2, for 
baselines between 25 and 200 m. A convergence 
test to the limit H — > oo is indicated for b — 200 
m. 



the remote, high-flying layers of air. As we turn 
to the astrometric task of completing the right- 
angled triangle formed by D, P and b* in Fig. 3 — 
with the aim of precise determination of either of 
the base angles at b* — , the "baseline calibration" 
emerges as an additional focus. This means deter- 
mination of b* , the image of b. The following sec- 
tions consider the atmospheric lensing correction 
b* — b to a, in principle, rock-solid and accessible 
ground baseline b. The computational strategy is 
to derive P and D, then to use Eq. (7). Within 
this framework, the geometric baseline b is defined 
joining the "ends of the paths through the atmo- 
sphere," and is assumed to be the same vector for 
both stars in the case of astrometry. We do not 
ask the question whether the corresponding ter- 
minal points Ti and T 2 are that well defined for 
real telescope optics with chromatic, azimuth de- 
pendent foci and trusses that bent under the load 
of their own weight or the wind. 

The formal solution to the problem of tracing 
four beams (two stars, separated by an angle r and 
labeled P for "primary" and S for "secondary," to 
two telescopes) for a given Z in Eq. (6) is then 
given by computation of Dp for the primary at 
some Z (some baseline b, see Eq. (5)), computation 
of Pg for the secondary at the same Az = Z = 
z iS ~ Z 2S = z iP ~ Z 2P ( smce * ne secondary must 
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be caught by the same two telescopes) but slightly 
different true zenith angles z 1 g = z 1 p + r, z 2 g = 
z 2 p + t. It should be noted that the earth-bound 
baseline is not tilted here — there is no lifting of 
one telescope and sinking of the other to acquire 
the secondary — , and b refers to the geometrical 
distance between two foci that define a common 
reference for both stars. Any residual effects of the 
spherical atmosphere and/or atmospheric layering 
are then caused by the second derivative of L(z) 
(Fig. 8). 




Az 



Az 



Fig. 8. — A difference in the delays Dp and Dg 
measured for the primary star and the secondary 
star is caused by the nonlincarity of the path 
length L, Eq. (14), as a function of the true zenith 
position. 

From (23) with (22), omitting a constant p + H, 
L = - J (p + H) 2 - p 2 n 2 sm 2 z 



P+H 



r—p 



(r 2 n 2 



TTi 

2 2-2 U / 2 <M 26 ) 

P z Hq sin zq) 1 '^ 



where we insert zq = z — R. For a spherical earth 
without atmosphere, the vacuum limit, 

R n -=l0, no— >1, (27) 

L — ► p(l — cos z) 

p 3 5 

= -tan 2 z p tan 4 zH ptan 6 z.(28) 

2 8 16 v ; 



D — ► p(cos z 2 — cos z\). 



(29) 



This suggests we define an effective true zenith 
angle z via 



D = bsinz — Ptanz, 



(30) 



and use Eqs. (5) and (29) to prove that this equals 
the mean, 

(3D 



z = 



Since the astrometric signature is hidden in the 
second derivative of L(z), the computationally 
most appealing expansion of R and L is a Tay- 
lor series around some reference value of z: 

R\z+x = £,o+£,ix+&x 2 +&x 3 + . . . , £ = R\ z - 

(32) 

We define refractivity integrals as a short-cut to 
the notation, covering Eq. (8) as a special case: 



R 3 



= /2j+l 
/ = 



no 



dn 



(33) 



n(r 2 n 2 - / 2 )J+i/2 : 
pn sinz , j = 0,1,2,... (34) 



A stable numerical scheme for these integrals is 
proposed in App. B. Insertion of the series (32) 
into the l.h.s. of (8) and into the arguments z = 
z — R of the sines at the r.h.s. yields the expansion 
coefficients 



Co — Ra, 

Ro + R\ 



6 



tanzo 



/ 1 + 



i?0 + R\ 

tanzo 



Rq + i?i — 3 



R\ + i?2 



tan zo 



[3i?i2 — Rqi tan 2 zq] tanzo 



2i& 



(35) 
(36) 



(37) 



-(6 - I) 5 



(i?o + i?i)6 + 3^^6 



+ 



1 a H ° 

;an zn Lb 



tan zo 
1 i?i 



- 3- 



tan 2 zo 

R2 

tan 2 zo 



5 Ra 



2 tan 2 z 2 tan 2 z 
tanZo ^5(^-^3)^01 



+ 



§Jfc}] (38) 



+9i?i 2 (3i?i2 



2P i - -Roi tan 2 z ) 



+i?oi (-Roi + 3i?oi tan 2 z ) tan 2 z , (39) 
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with the doubly- indexed shorthands 

R l3 = R l + R J , (40) 
Rij = Ri + R 3 ■+ tanz , (41) 
Rij = Ri + Rj ■— tanz . (42) 

In Eqs. (33)-(42) and App. B, the subscripts of R 
are the exponential j of the definition (33); else- 
where they indicate the telescope number/site, x 



0.0020 I 1 1 1 1 1 r-| 

0.0018 - / - 
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0.0002 - Z'^^^L^^ " — ° 

0.0000 I — : — — t~~ < — 1 1 1 
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true zenith angle z (deg) 

Fig. 9. — The expansion coefficients £j of Eq. (32) 
as a function of zenith angle z for j = 0, . . . , 3 and 
the atmospheric model (11)— (12). R is an odd 
function of z, so the even- indexed 6 approach 
for z — > 0. £i approaches n — 1, see Eq. (1). 

in Eq. (32) is of the order of b/(2p) if the ref- 
erence azimuth z is chosen close to the middle 
between the telescopes, and therefore not larger 
than 1.6 • 10~ 5 rad for b < 200 m. Because the 6 
are approximately of the same magnitude (Fig. 9), 
collecting the terms up to j = 3 ought establish 
a relative accuracy of w 5 • 10~ 14 in the angle of 
refraction. 

The expansion 



■>\z+x 



+ E 

i=l,2,3,. 



(43) 



proceeds via insertion of Eq. (32) into the sines of 
Eq. (26), and employs an auxiliary set of integrals 



[{p+Hf-py- 1 / 2 

rP+H 
+(2i-l) / 

J r—p 



(r 2 n 2 - J2)*+i/2 



dr 



— v p tan 1 " z cos z . (44) 



h 
h 



1 £l n— >psinz 

v\ — ► psinz, 

tan zq 



-2 Vl 



(1 



1 )d-6) 2 + 26 



tan^ zq 



tanz 



i (1-6) 2 

-V2- 



2 tan^ z 

n—>l,I—>-psinz 1 
1 



^l(l-^l) 

1 1-6 
-— 1>2 

2 tan zq 



— pcosz, 
■2 ^—^ +3g 2 (l 



(45) 

(46) 
(47) 



(1 



tanzo 
1 



5—) 

tan zo 

26 



----- )(i-6) i 

tan z z tan zo _ 



1^ (1-gl) 3 n-1,7- 

2 tan 3 zq 



• psin 3 



l 



7/9 sin z, 



ra— >1,J— *psin z 1 



pcosz. (48) 

24 P v ; 

If the li are calculated at the mean position 
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Fig. 10. — The differences between the Taylor ex- 
pansion coefficients lj of Eq. (43) and their vac- 
uum values as a function of zenith angle z for 
j = 1, . . . , 3, parametrized through the exponen- 
tial model (11)- (12). L is an even function of 



z, so the odd-indexed L approach for z 



0. 



The vacuum limits, Eq. (29), are h 
h — * /o(cosz)/2, l 3 — ► — p(sinz)/6. 



psmz, 



z = (zi + z 2 )/2, the delay for a single sky posi- 
tion becomes 



D w h (zi - z 2 ) + l 2 (z\ - 4) 

+h {-4 - zi) + u (zf - 4) + ■ ■ ■ (49) 

w / 1 Az + ? 3 (Az) 3 /4 + ? 5 (Az) 5 /16 

+Z 7 (Az) 7 /64+--- (50) 
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If the U are calculated at the mean position z — 
(z 1 p + z 2 p)/2, the differential delay is to lowest 
orders in r 



D s -D p 



Azr 



2Z 2 + 3rZ 3 + ((Az) 2 



4r 2 ) h 



+ (^(Az) 2 t + 5t 3 J l 5 

+ Q(A^) 4 + 5(Az) 2 r 2 + 6r 4 ^) l 6 + ■ ■ -(]51) 

The leading term 2Azrh contains 

• the familiar geometric "vacuum" contribu- 
tion Azprcosz 6tcosz; see Eqs. (6) and 
(47), 

• an atmospheric correction to fo of the order 
of 1-10 m/rad 2 (Fig. 10). It adds some tens 
of nanometers to the differential delay, if r < 
1 arcmin = 3 • 10~ 4 rad, and Az < 3 • 10~ 5 
rad (b < 200 m) . There is no equivalent con- 
tribution of this kind in planar earth models 
like Fig. 1. 



h 



With Eqs. (15) and (32), 
prig sin(zi — R\) 
pn a sin (z- £ + (1 - £i] 



Az 



(Az) 2 



(Az) 3 



pn a sin(z 2 - R 2 ) 

pn Q sin (z - Co + (£i - 1) 



Az 



-i 



(52) 

). (53) 
(54) 
(Az) 2 



(55) 



we may expand P 2 in a power series of Az, 
P 2 « {pn cosz (l- a)} 2 (Az) 2 

x {[(6 - l) 3 - 66] coszo + 66 [1 - £i] sinz } 
x(Az) 4 + ... (56) 

We do the same for D 2 via Eq. (50), and eventually 
combine these power series of Az at the r.h.s. of 
Eq. (7), 

b* 2 ^ {(I -^) 2 p 2 n 2 cos 2 z + l 2 ){Az) 2 



+ 



f p 2 nl(l-^ 1 )cosz / 3 

| ^2 V^ 1 ~ ^ ^ ° OS Z ° 

+66[l-ei]sinz )+^|(Az) 4 



(57) 



which turns into Eq. (5) in the vacuum limit. Ex- 
amples of this "baseline magnification" introduced 
by the atmosphere are shown in Fig. 11; the effect 
becomes larger if the transition into the free space 
is smoothed by choosing a large cut-off height H. 




20 30 40 

z, (cleg) 

Fig. 11. — The difference b* —b between the length 
of the effective baseline above the atmosphere and 
the geometric baseline on the ground according to 
Eq. (57) as a function of z\, parametrized through 
the model (11)-(12), which is cut off at H = 20, 
40 or 80 km. 



2.2. 3D Geometry 

2.2.1. Geographic Coordinates 

The previous section dealt with the case where 
the telescopes, the star and the earth center are 
coplanar. In the general case, the direction of the 
star and the direction of the second telescope do 
not share the same azimuth A. Let the telescopes 
have geographical latitudes <&i and longitudes Aj 
in a geocentric spherical coordinate system: 



n = p 



cos A, cos $i \ 

sinAjCos4>i , (i = 1,2). (58) 
sin $j / 



For the Very Large Telescope Interferometer in 
Northern Chile $ w -0.4298 rad and A w -1.229 
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rad, for instance. The baseline angle in Eq. (3) 
becomes 

cos Z — sin $1 sin $ 2 + cos $1 cos $2 cos(AA), 

(59) 

where A A = Ai — A 2 . To transform Cartesian 
coordinates from the local alt-az-system of tele- 
scope i (with the Cartesian coordinate z pointing 
to the zenith, x horizontally tangentially to the 
earth toward north and the local horizon as in- 
dicated by the dotted line in Fig. 2, y horizon- 
tally toward west) to the geocentric system (with 
z pointing from the earth center to the north pole, 
x from the center to the equator south of Green- 
wich, y from the center to the equator 1000 km 
west of Sumatra) we translate the coordinates into 
a tilted system originating from the earth center, 
then (de)rotate them: 



— sin <f>i cos Ai sin Aj 

— sin $j sin Ai — cos A., 



cos $ 







cos $i cos Ai 
cos $i sin Ai 
sin <I>i 




i = 1,2. 



(60) 



The inverse operation with the inverse matrix 
(which equals the transpose matrix) is 

— sin $i cos Ai — sin $i sin Ai cos $i 

sin Ai — cos Ai 

cos $i cos Ai cos <3>i sin Ai sin $i 




(61) 



The product of two such operations with indexes 1 
and 2 converts the two alt-az systems: starting in 
Eq. (60) with r 2 = 0, computing r c , and inserting 
this into (61) for i = 1 shows that the origin of 
coordinates of telescope 2 is located at 

bi2 = p ■ 

— sin $1 cos $2 cos(AA) + cos <3>i sin $ 2 

cos $ 2 sin(AA) 
cos $1 cos $2 cos(AA) + sin $1 sin $2 — 1 

seen from the origin of telescope 1. The length of 
this vector is b = |bi 2 | of Eq. (3), using Eq. (59); 
the third coordinate is negative since the second 
telescope lies below the horizon of the first tele- 
scope (and vice versa), as illustrated in Fig. 2. 



2.2.2. Vacuum limit 

If the atmosphere is absent, the star direction 
is defined as 



cos Ai sin Zi 
sin Ai sin z% 

COS Zi 



1,2, 



(62) 



in terms of the true local azimuth Ai and true 
zenith angle zi in the zth telescope coordinate sys- 
tem. A is counted positive starting from N to 
W. [This azimuth convention is the one of (Smart 
1949, §11); the alternative convention of (Taff 1980; 
Karttunen ct al. 1987) is obtained with the re- 
placement Ai — > it — Ai] The cosine of the angle 
between the star and the baseline in Fig. 2 is 



cost^i = si • bi 2 /6, 



(63) 



where Si • bi 2 is known as the geometric optical 
path delay. If one swaps the indexes 1 and 2, the 
cosine switches its sign, because in this parallax- 
free situation the angle between star and baseline 
is the 180°-complemcnt of the angle relative to the 
other telescope: 



cos</?2 = S2 • h 2 i/b = — cos<pi. 



(64) 



This may be verified with the standard coordinate 
transformations between the hour angles hi and 
right ascension 5 for i = 1,2, 



cos S sin hi = sin Zi sin Ai , 



(65) 



sin Zi cos Ai = cos $i sin S — sin $i cos S cos h i7 

(66) 

cos Zi — sin $i sin S + cos $i cos <5 cos hi , (67) 



hi - h-2 = Ai - A 2 



(68) 



So if the atmosphere is absent, this angle relates 
to P = I sin (fi\b (i = 1,2) as shown in Fig. 3. 

The mean and difference in the true zenith an- 
gles remain defined as in Eqs. (31) and 



Az = zi — z 2 , 



(69) 



and can be retrieved from the geographical coordi- 
nates (66)-(67) and (Abramowitz & Stegun 1972, 
(4.3.34)-(4.3.37)). 
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2.2.3. Ray tracing 

There is one distinguished geocentric coordi- 
nate system for the case of a single star, shown 
in Fig. 12, in which the direction vectors of the 
incoming light above the atmosphere have only a 
component along the polar axis. The side view 
of this geometry reduces to Fig. 6 if the posi- 
tions 1, lv, 2 and 2v lie on the same projected 
straight line. (In the following, "projected" means 
projected onto a plane perpendicular to the ray 
propagation above the atmosphere.) The strat- 




Fig. 12. — A top view of the earth coordinates as 
seen from the incoming two rays in which both 
Si are perpendicular to the plane of the draw- 
ing. If the atmosphere were absent, the first ray 
would see telescope 1 at the position lv. The at- 
mospheric refraction pulls rays toward the earth 
center, which means it has to relocate its impact 
at the top of the atmosphere to the actual image 
position 1 above the atmosphere to end up at tele- 
scope 1 on the ground. The same effect for ray 2 
means the projected baseline vector P above the 
atmosphere from 1 to 2 is both longer and tilted 
by Aij compared to the vacuum case from lv to 
2v. 

egy to transform the star position to the projected 
baseline (vector) P is (i) to calculate the trian- 
gle formed by the vacuum baseline and the earth 



center seen by the incoming light, and (ii) stretch 
this radially outward to include the effect of the 
two difference refraction angles Re. the two true 
zenith angles Zi for both telescopes are assumed 
to be given via Eq. (67). The impact parameters 
Ii of the rays for the vacuum case are p sin Zi and 
they equal the lengths of the (projected) station 
vectors from the earth center to the points lv and 
2v in Fig. 12. The vacuum baseline in Fig. 12 
stretches from lv to 2v, which has the projected 
length b\ sin<^| as written down in Eq. (63). The 
projected baseline aperture angle a in the trian- 
gle with side lengths p sin Zi and 6|siny>j| relates 
to the impact parameters and projected vacuum 
baseline by planar trigonometry (Abramowitz & 
Stegun 1972, 4.3.148), 



sin 2 z\ + sin 2 z 2 



2 sin z\ sin zi cos a H — ~ sin 2 w\ . 

P 2 

(70) 

a is the projection of Z defined in Eq. (59). The 
equation mingles the true zenith angles from the 
two telescope's point of view with the angles ipi, 
which represent the star distance from the baseline 
direction. One may reduce this to the geographic 
coordinates by multiplying ri + bi2 = r 2 with s\ 
to get 



coser 



COS Z — COS Z\ COS Z2 



= 1 - 



sin z\ sin z-i 

COs(£i - Z 2 ) 



cos Z 



sin zi sin Z2 



(71) 
(72) 



cos a and cos Z are close to unity in contemporary 
optical interferometry, so the actual implementa- 
tion avoids the use of the cosine in favor of the 
haversine, 



haver 
Z 2 



sin z+^z sin z^z 



sin z\ sin Z2 

16 48 



(73) 



4 sin 2 z 



■ sin 2 z 



Az 



smz 



1 



48 Vsinz 



1 Z 2 {Azf 
16 sin 4 z 



(45 - 30 sin 2 z + 2 sin 4 z) 



{Azf 



+ (3 - sin 2 z) 



+ - 



Z & 



1440 sin 2 z 



2880 sin 6 z 
{AzfZ 2 {Az) 2 Z A 

192 sin 6 z 192 sin 4 z 
+ ..., 



(74) 
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to protect against cancellation of significant digits, 
and returns from there to the sine, if needed, 



sin cr = 2-v/havoVl — haver (75) 
(2j-3)!! /havcr XJ 



2v / hava- 



Vhav cr ( 2 — hav cr — 



hav cr hav cr 



The special coplanar case of Sec. 2.1 is included as 
Az = Z and a = or tt. (On some Linux systems, 
where the cosl library function does not sup- 
port the full accuracy of the 96 bit long double 
number representation, it makes sense to switch 
to alternative high-precision implementations of 
the cosine (Schonfelder 1980).) Two calculations 
for the actual impact parameters 1^ — pn sin zjp 
starting from the given Zj_ would be done as in 
the preceding, "aligned" geometry of Sec. 2.1, and 
these inserted into 



Il+Il = 2/i/2 COS CT + P 2 



(76) 



to calculate the projected baseline P and to gen- 



eralize Eq. (15). P 2 = (h - h 



4/1/2 haver 



comprises the terms of Eq. (56) of the constrained 
geometry augmented by 

4/1/2 hav a = (pn Q sin z ) + ... (77) 

sin z 

as read from Eq. (74) in combination with Eqs. 
(53)-(55). 

To calculate the path difference D = Li—L 2 be- 
tween the two rays above the atmosphere, Eq. (18) 
remains valid, but in general a Taylor expansion 
akin to (19) does not exist, because P > li sin a 
(i = 1,2) excludes small P at arbitrary a. Eq. 
(32) remains in use to expand the refraction angle 
in a neighborhood of z, and so do Eqs. (43)-(50) 
that depend only on zenith distances but not on 
azimuths. 

2.2.4-. Baseline Rotation 

The rotation angle between the baseline b (pro- 
jected on a plane perpendicular to the star direc- 
tion) and P is A77 = r/p — rjb where rjb and rjp are 
the angles from b to / 2 and P to / 2 respectively, 
and given by 

h 



h 



= cos cr + sin cr cot rjp, 



(78) 



sin z 2 
sin Z\ 



= cos a + sin tr cot 



(79) 



Numerical examples are presented in Figs. 13- 
14: By symmetry, the rotation effect vanishes if 
the star azimuth is along the baseline or perpen- 
dicular to it. The angle A77 may be about five 
times larger than the intcrferomctric resolution of 
a 200 m baseline in the K-band: the interferomet- 
ric fringes appear slightly rotated on the detector, 
and the true (u, v) coordinate is found by rotating 
the "apparent" vector by — A77. 




100 110 120 130 140 150 160 
P(m) 



170 180 190 200 



Fig. 13. — The baseline rotation A77 introduced 
in Fig. 12 for a baseline b = 200m, for 6 differ- 
ent azimuth angles Ai measured from T\ toward 
the baseline, and for sinzi changing in equidistant 
steps of 0.05 from 0.05 to 0.85. 

There is no intrinsically new aspect compared 
to the analysis of Sec. 2.1: The rotation of the 
baseline vector P relative to the projection of b 
would again be absorbed into the pointing direc- 
tion of the two telescopes. With only one (scalar) 
observable, which is the differential delay, one can- 
not measure both the astrometric angle (distance) 
between two stars and the positioning angle of the 
S relative to the P at the same time. 

2.2.5. Astrometric Case 

Eq. (51) obviously looses its meaning because 
the secondary star now has two degrees of freedom 
and can no longer be positioned by a single angle 
t. The details follow by writing down Eq. (67) 
for both telescopes and both stars. If the star 
distances th and t$ are defined as 



h s = h P + T h , 5 S = dp + t s , 



(80) 
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z, (deg) 

Fig. 14. — An alternative view on the six lines of 
Fig. 13 with the abscissa switched from the pro- 
jected baseline length to the azimuth angle. Each 
small cross has a counterpart in Fig. 13. This 
graph would not change visibly choosing a baseline 
of b = 100 m. 



with numerical and analytical methods known 
from treatments of the more familiar refractivity 
integral. Subtraction of two of these computes the 
optical path difference, and renormalizes the op- 
tical path lengths to start from a common plane 
tangential to the earth's upper atmosphere, which 
at the same time defines the projected baseline in 
this plane perpendicular to the two rays. Defini- 
tion of a right triangle above the atmosphere with 
the projected baseline and the path difference as 
two catheti defines a hypotenuse, which is an ef- 
fective baseline that is longer than and rotated 
relative to the geometric baseline between the re- 
ceiving telescopes on earth. 



the expansion of Eq. (67) yields for z.^ = z-ip + Ti Z 
to lowest order in tk and Tg: 

cos $i sin 5p cos hp — sin <&j cos <5p 

Tiz = : T$ 

sin Zip 

cos $i cos <5p sin hp , „ . 

H : r h + ... (81) 

sin zip 

The changes tm in the azimuths are not detailed 
here, since the refraction is determined by the 
zenith angles; they are expected to ensure that the 
associated change in the star direction (62) forms 
an acute angle to the baseline to maintain the in- 
terferometric resolution. The symmetry suggested 
in Fig. 8 will generally be broken: Azp ^ ^zg. 
Restarting from Eq. (50), the differential delay is 
expanded in powers of the doubly differential At z : 



= Az P f z 



2l 2 + 3fJ 3 + ((Azp) 2 + 4f z 2 ) k 



+ ■ 



+ 0(At z ), (82) 



with t z = (t 1z +t 2 2)/2 and Ar z = t 1z -t 2z , where 
the lj are again evaluated at the mean primary 
zenith, (z iP + z 2 p)/2. 

3. Summary 

The optical path length integral of star light 
passing through the atmosphere can be handled 
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A. Atmospheric Refraction Models 

For the model (11), the misassignment in the star position is a few mas for each kilometer of error in the 
scale height K. In the limit K —> 0, we recover the values for the flat earth of Eq. (1); the corresponding 
mathematical argumentation is given by (Stone 1996, (4)) and the limit H* — > in (Livcngood et al. 1999, 
(1)). In practice, the scale height is coupled to the atmospheric gas density, and in a self-consistent model of 
a single average molecule species like Eq. (11), K is uniquely coupled to the atmospheric pressure at ground 
level (McCartney 1976). 

In a variant of this problem, the scale heights of various gas components differ and their mixing ratios 
change. In this case, Eq. (11) would be upgraded to a sum over gas components with individual pairs of Xo 
and K — it is doubtful to scale the entire, single scale height with the ground humidity instead (Livengood 
et al. 1999). Consider for instance a fixed, "dry air" contribution to %(r) with a constant scale height of 10 
km at xo = 4.0808- 10~ 4 , superimposed by water vapor at xo = 1-75 • 10~ 6 that "freezes out" at scale heights 
of only a few kilometers: the pointing variations are w 25/xas per km change in the water vapor scale height. 
Calibration of this water column could be achieved through monitoring spectral ranges of high atmospheric 
water absorption (Akeson et al. 2000; Meisner & Le Poole 2003). 

The small difference in the refractive indexes at wavelengths of 2 and 2.4 /im, about A%o ~ 2.1 • 10~ 7 , 
displaces these two colored rays horizontally by about 0.80 mm if they hit the earth surface, calculated at 
apparent zenith angles of zq ss 30 deg. If one looks at the starlight as an unvignetted plane wave, there is no 
such effect. Instead, there is a rainbow effect caused by this lensing of the earth atmosphere, which spreads 
the apparent positions of these two colors on the sky by 12.5 mas, which can be estimated through Eq. (1) 
as AR w tanzoAn using An Axa/2 with Eq. (10). The calculation within the ray optics, however, can be 
used to consider the distortion/decorrelation induced by turbulence on length scales of these displacements 
(Colavita et al. 1987). 



B. A Numerical Approach to the Refractivity Integrals 



In the course of this investigation, the integrals (33) have been decomposed into refraction within the 
open interval p < r < p + H, plus one kink at r = p + H where e = n 2 changes abruptly from 1 to e r=p+ n ■ 



pj+i 
2 

pj+i 



r 



de 



e(r 2 e - 7 2 )J+ 1 /2 
de 

e(r 2 e _ /2)J+ i/2 



+ I 2 ^ 



n r=p +H 



dn 



i[(p + H) 



py+i/2 ■ 



(Bl) 
(B2) 



The second term is Snell's law in terms of the two angles of incidence, simpn = I/(p + H) and 
n r=p+ H sin^fl_ = sinipH just above and below the atmosphere top boundary, 



n r=l ,+H 



dn 



( P +g) 2 » 2 



- l)j + l/2 



JlpH- 



tan 2j xdx = < 



IpH - IpH- 

tan(^jy) -i> H - 



,3=0 
,3 = 1 



tan (ip H ) 



1 



tan(V>jr) + ipH - ■ ■ ■ ,j = 2, 



(B3) 

where the three dots mean the previous expression is to be subtracted with all ipn replaced by ipH-- The 
first term in Eq. (B2) is mapped onto the interval < t < 1 through the substitution 



i = 7 



r — p 
r - p + H(l-~f)' 



(B4) 
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7 is chosen with the idea of important sampling such that r\ /2 — typically selected as 2 km — is mapped on 
t = 1/2. This yields 7 = (H - r 1/2 )/(H - 2r 1/2 ) and leaves 



Je r=p+H t{r 2 t - PV+ 1 ' 2 = ~ X) A (i - 7) 2 e(r2e - /2)j+i/2 ' ( B5 ) 

This has been integrated with the trapezoidal rule, doubling the number N of equidistant sampling points 
after each step. As e(r 2 e — I 2 y +1 / 2 is a smooth function over the t-interval, the error in the trapezoidal 
rule is dominated by the variation in (dx/dr)/(t — -f) 2 . For the exponential model (11) one can demonstrate 
that the sequence Vat, Vzn, Van,- ■ ■ , obtained by repeated division of the sampling step size by two, obeys 
Van — V2N ~ |(^2JV — Vjv). The Richardson extrapolation induces the estimator V N Zl?° 2V2N — Viv, not 
the Simpson rule; monitoring the convergence of this estimator has been used to terminate the subdivision 
loop. A very similar analysis is applicable to the integrals Uj. 
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